# Get Table 5
# =========================================================================================================
# remove everything
rm(list=ls())

# =================================================================================

# Load packages
library(foreign)

# Load the bootstrapped data
load("Boot_s_Dt.RData")
load("Jack_s_Dt.RData")

# ---------------------------------------------------------------------
# Get the bootstrapped CI.
B <- 400
dec <- 4

Bcoef[c(4,5,8),1]; Bcoef[c(7,8,9),1]

# acceleration factors
acc.top <- apply((apply(JKcoef, MARGIN = 1, mean) - JKcoef)^3, MARGIN = 1, sum)
acc.bot <- 6*(apply((apply(JKcoef, MARGIN = 1, mean) - JKcoef)^2, MARGIN = 1, sum))^(3/2)
acc.factor <- acc.top/acc.bot

# CI
for (pos in c(4,5,8)) {
  phio <- qnorm(sum(Bcoef[pos,] < Bcoef[pos,1])/B)
  phia <- qnorm(0.025); phina <- qnorm(0.975)
  # a_1 and a_2
  c <- acc.factor[pos]
  a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
  a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
  up <- format(round(quantile(Bcoef[pos,], probs = c(a1)), dec), nsmall = dec)
  down <- format(round(quantile(Bcoef[pos,], probs = c(a2)), dec), nsmall=dec)
  ci <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
  print(ci)
}

for (pos in c(7,8,9)) {
  phio <- qnorm(sum(Bcoef[pos,] < Bcoef[pos,1])/B)
  phia <- qnorm(0.025); phina <- qnorm(0.975)
  # a_1 and a_2
  c <- acc.factor[pos]
  a1 <- pnorm(phio + (phio+phia)/(1-c*(phio+phia)))
  a2 <- pnorm(phio + (phio+phina)/(1-c*(phio+phina)))
  up <- format(round(quantile(Bcoef[pos,], probs = c(a1)), dec), nsmall = dec)
  down <- format(round(quantile(Bcoef[pos,], probs = c(a2)), dec), nsmall=dec)
  ci <- paste(paste("(", up, sep = ""), paste(down, ")", sep = ""), sep = ", ")
  print(ci)
}

